Universal energy fluctuations in thermally isolated driven systems 
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When an isolated system is brought in contact with a heat bath its final energy is random and 
follows the Gibbs distribution - a cornerstone of statistical physics. The system's energy can also be 
changed by performing non-adiabatic work using a cyclic process. Almost nothing is known about 
the resulting energy distribution in this setup, which is especially relevant to recent experimental 
progress in cold atoms, ions traps, superconducting qubits and other systems. Here we show that 
when the non-adiabatic process comprises of many repeated cyclic processes the resulting energy 
distribution is universal and different from the Gibbs ensemble. We predict the existence of two 
qualitatively different regimes with a continuous second order like transition between them. We 
illustrate our approach performing explicit calculations for both interacting and non-interacting 
systems. 



Understanding equilibrium and non-equilibrium prop- 
erties of thermally isolated systems has become a fore- 
front of research due to experimental developments over 
the past decade, particularly in cold atom systems [J, 
trapped ions Q , and nuclear spins 0] and superconduct- 
ing qubits In these systems the coupling to external 
dissipative degrees of freedom is strongly suppressed and 
irrelevant on accessible time scales. These systems pro- 
vide a new and very clean playground where one can in- 
vestigate fundamental questions in statistical and quan- 
tum physics. Moreover, they point to new practical ap- 
plications, in particular in the context of quantum in- 
formation. The experimental studies inspired intensive 
theoretical research on a variety of topics. These include 
equilibration in isolated systems initially driven out of 
equilibrium by a sudden change in a coupling constant (a 
quench) ; defect (or energy) generation during slow nearly 
adiabatic processes in gapless phases or near singulari- 
ties, such as quantum phase transitions (for a review see 
Rcfs. Hi); non-equilibrium quantum phase transitions 
in the presence of 1// noise [7|; and many more. 

In this work, we consider the energy distribution of a 
thermally isolated system following a non-adiabatic pro- 
cess. Consider two setups where the energy of an isolated 
system is changed. In the first, the system is brought into 
contact with a heat bath until equilibration and then dis- 
connected from it - similar to an oven. In the second the 
energy of the system is increased due a non-adiabatic 
change of some external parameter(s) - much like a mi- 
crowave. The two setups are illustrated schematically in 
Fig- [II Our interest is in the energy distribution of each of 
the systems at the end of the process. For the first setup 
the result is well-known and corresponds to the classic 
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heating mechanism which can be found in any book on 
thermodynamics (see e.g. Ref. Q). If, as usual, the bath 
is large compared to the system then the energy distri- 
bution of the system becomes Gaussian, with a canonical 
width uniquely determined by the fluctuation-dissipation 
relation; SE"^ — T^Cy, where T is the temperature and 
Cy is the specific heat. This relation is valid for both 
quantum and classical systems and is independent of the 
details of the interactions between the system and the 
bath. Now, consider a second setup where the energy 
of the system is changed due to a non-adiabatic varia- 
tion of an external parameter (say, the electro-magnetic 
field in the case of a microwave or the motion of the 
piston in Fig. [1]). This type of heating also inevitably 
leads to an uncertainty in the final energy of the system. 
While it is known that the energy changes in the system 
obey, even beyond linear response, the recently discov- 
ered fluctuation theorems (see for example, j9l-[ll|) very 
little is known about the resulting energy distribution 
in this case [T^]- In a large macroscopic system the en- 
ergy distribution is expected to be very narrow and the 
relative energy fluctuations negligible. But in small or 
mesoscopic systems, which are of primary experimental 
interest (see e.g. Ref. [ij), the fluctuations can be large 
and important. Fundamental questions are unanswered: 
Which features of the energy distribution are universal 
and which features depend on details of the system and 
driving protocol? To what extent can the width of the 
distribution be controlled? For example, can one dynam- 
ically increase the energy of an isolated system without 
increasing the uncertainty in the final energy? Can the 
fluctuation-dissipation relations, which determine the en- 
ergy width in the oven-like setup, be extended to the 
microwave-like setup? 

In this paper we begin to address such questions. 
Specifically, we study a thermally isolated system un- 
dergoing a repeated cyclic process, whereby some exter- 
nal parameter A in the Hamiltonian is changed in time 
and returns to its original value at the end of each cycle, 
see Fig. [T] We show that under generic assumptions of 
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FIG. 1. Schematic comparison between the usual thermal 
heating (traditional oven, top) and an energy increase due to 
not-adiabatic work (microwave oven, bottom). On the right 
of each case we present a schematic picture of the resulting 
energy distribution in each case. 



(i) small work per cycle and (ii) absence of correlations 
between cycles (see detailed discussion in Methods) the 
variance of the energy distribution (E) at energy E as- 
sumes a particularly simple form to leading order in l/N^ 
where N is the number of degrees of freedom in the sys- 
tem. It depends only on the microcanonical temperature 
/?(£') = dE^nil{E) where n{E) is the density of states, 
and on the average energy change in a cycle at energy E, 
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Here Eq is the initial energy of the system and CTq is the 
initial variance. This equation is the main result of the 
paper. 

The result ^ follows from integrating a Fokker-Planck 
equation which describs the time evolution of the en- 
ergy distribution P{E^ t) (see Methods and Supplemen- 
tary Material for details): 



1 



dtP = -dE{A{E)P) + -dEE{B{E)P) 



(2) 



The change of the energy distribution in one cycle of the 
protocol is obtained by integrating this equation over 
the duration of the protocol, set for simplicity to be 
unity. Within this choice A{E)^ B{E) represent the aver- 
age work per cycle and its variance respectively: A = {w) 
and B = {w'^)c- Here the angular brackets denote aver- 
aging over realizations of the cycle starting from a fixed 
initial energy. 

In general A{E) and B{E) are protocol dependent 
functions and are a priori independent from each other. 
However, since the system is thermally isolated its time 
evolution is governed by Hamilton's equations of motion 
in the classical case and the Schrodingcr equation in the 
quantum case. This puts strong constraints on the re- 
lation between A{E) and B{E) similar to the Einstein 



fluctuation-dissipation relations between drift coefficient 
(mobility) and diffusion in open systems [l3|. In partic- 
ular, we find 



I3B^2A- OeB = 2A + 0(A^"^). 



(3) 



For interacting systems with many degrees of freedom the 
second term on the RHS of this equation is a 1 /N correc- 
tion which can be neglected. This term can be important 
though in mesoscopic or integrable systems. Eq. ([3]) was 
previously suggested for classical systems in Ref. [2^ . 
The fluctuation-dissipation relation ([3]) is derived within 
a small work assumption (explicitly (3{E)'^{w'^)c <C A{E), 
where {w^)c is the third cumulant of the work, and 
A{E) <C TCy, see Methods and Supplementary Mate- 
rial) . As we will show below Eq. ([3]) holds for a very wide 
class of classical and quantum systems starting from non- 
interacting particles in a time dependent cavity to fully 
interacting spin systems. The main result of the paper 
Eq. ([T]) is a direct consequence of this relation (see Meth- 
ods). 

Several interesting consequences follow from Eq. ([Ij: 
(i) When A{E) is constant the energy width depends 
only on P{E), and not on the amplitude of the drive or 
other details of the driving protocol, (ii) When A{E) is 
not constant, depending on the functional form of A{E) 
and 13(E), the variance of the distribution can be larger 
and surprisingly, even smaller than the width of the equi- 
librium Gibbs distribution at the same mean energy. In 
fact, a"^ [E) / (j1^{E) can be made arbitrarily small by a 
proper choice of A{E). (iii) When A is a function of 
the energy density u ~ E/N (with a possible extensive 
energy independent prefactor like the total number of 
particles), we have cr'^{E) ^ 0{N), scaling as in equilib- 
rium. For a single quench this result was noticed e.g. in 
Ref. • Here we show that it remains valid after many 
quenches, (iv) The dependence of on E displays two 
qualitatively distinct behaviors with increasing E, de- 
pending on whether the integral in Eq. ([1]) diverges or 
converges as E ^ oo. 

To illustrate the distinct behaviors associated with 
point (iv) above we consider the generic case where 
/3 cx E~", which is the case for phonons, superfluids or 
other systems with Goldstone bosons, Fermi liquids, ideal 
gases and others. Moreover we measure time in units of 
the number of cycles carried out (in what follows we will 
use time and number of cycles interchangeably) and as- 
sume a simple power law behavior for A{E): 



dtE = A{E) = cE'' 



(4) 



As will become clear below the two regimes exist even 
in cases when A{E) and P{E) are not power laws. The 
values of a are constrained by simple thermodynamic ar- 
guments to < a < 1: the lower bound is required by 
positivity of the specific heat and the upper bound as- 
sures that the entropy {S{E) cx E'^"") is an increasing 
unbounded function of the energy. To prevent the sys- 
tem's energy from diverging in a finite time we require 
s < 1 (as follows from integrating Eq. 
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For simplicity we also assume (Jo{Eq) = and com- 
pare the width to the equilibrium canonical width a^^ = 
—df}E ^ E^^°'/a. In this case the system displays a 
transition between two behaviors as the functional form 
of A{E) is changed. This transition is continuous and is 
characterized by a diverging time-scale needed to reach 
the asymptotic regime. Specifically, depending on the 
sign of 77 = 2s — 1 — a, Eq. ([T]) implies: (i) When 
77 < the width is Gibbs-like with (J^/cFcq ^ 2a/|77|, 
i.e. the ratio a'^ /a'^^ asymptotically approaches a con- 
stant value that can be either larger or smaller than one. 
Note that smaller widths correspond to protocols with 
large and negative s, i.e. to protocols where A{E) is 
a strongly decreasing function of energy, (ii) A second 
run-away regime occurs when 77 > 0. Here the width in- 
creases with a higher power of energy than the canonical 
width: CT^/(Tg^ ^ £■''. The resulting distribution is sig- 
nificantly wider than the canonical one. Given the con- 
straint on the value of s, this regime can only be reached 
if a < 1, in particular this regime is unreachable for a 
classical ideal gas. The transition between the Gibbs- 
like and run- away regimes occurs when 77 = which im- 
plies cF^/crlq ^ 2a In (^^^ ■ Close to this transition when 

1 7; I <^ 1, there is a divergent time scale (or number of cy- 
cles) required to reach the asymptotic regime. This time 
scale can be obtained by combining Eq. ([T]) and Eq. (U), 
see TableHl The crossover from the Gibbs-like to the run- 
away regime is qualitatively similar to a continuous phase 
transition, with the diverging time scale being analogous 
to a divergent relaxation time (critical slowing down) in 
the equilibrium case. We summarize our results, close to 
the transition, for the above choices of (3{E) and A{E) 
in Table m 

The qualitative difference between the two regimes can 
also be understood in terms of the entropy of the distri- 
bution: S = — p„ lnp„, where p„ are the microscopic 
probabilities to occupy different energy levels. Convert- 
ing this sum into an integral over energies and expand- 
ing the resulting expression up to corrections it is 
straightforward to check [l^ that: 



§iE) - §eq(£;) = In 



acqiE) 



where $cq{E) = \\i{\/2^:ac(^{E)Q.{E)) is the equilibrium 
canonical entropy. It is easy to see that the correction to 
the equilibrium entropy is always negative except when 
the width of the energy distribution coincides with the 
canonical width. In the Gibbs-like regime this correction 
is a constant, while in the run-away regime it has an 
explicit energy dependence. 

Note that under the assumptions used to derive our 
main result, the integral in Eq. ([1} can be rewritten in 
terms of A{S) the average entropy change per unit of time 
(per-cycle), where S = \aQ.{E). Using A{E) = A{S)dsE 
the nature of the transition between the two regimes as- 
sumes an interesting physical interpretation. Specifically, 
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TABLE I. A summary of the results for A{E) = cE\ f3 
E'" with < a < 1, s < 1 and 77 = 2s - 1 - a. The 
width specifies the asymptotic value in units of the equilib- 
rium width at the same energy. The time scale specifies the 
characteristic "relaxation" time needed to reach the asymp- 
totic regime. 



the integral in Eq. ([T]) becomes / dS/A^{S) showing that 
the Gibbs-like and (run-away) regimes correspond to the 
entropy growing slower (faster) than time squared. This 
conclusion does not rely on any assumptions about the 
specific functional form of 13{E) and A{E). 



I. EXAMPLES 

First, we consider a system of non-interacting and 
weakly interacting particles in a deforming cavity. Then 
we analyze a single particle in a harmonic potential, 
which is a part of a larger system, and subject to a time- 
dependent external force. Two additional examples of a 
classical one-dimensional XY-modcl and a quantum one- 
dimensional transverse field Ising model will be discussed 
in the Supplementary Material. 

Single Particle in a deforming cavity: Let us first con- 
sider a very simple system - a single particle bouncing 
elastically in a cavity. When the cavity is stationary 
the energy of the particle is conserved. If the cavity is 
chaotic there are no other conserved quantities so that 
in the long-time limit the particle relaxes to a uniform 
position distribution and an isotropic momentum distri- 
bution. We consider a process where the system is re- 
peatedly driven by deforming the cavity. At the end of 
each cycle the cavity comes back to its original shape 
and the system is allowed to relax in the sense described 
above (see FigH]). In this setup the number of degrees of 
freedom N is given hy N = 2d where d is the dimension- 
ality of the system. In a single collision with the moving 
wall the particle's kinetic energy can either increase or 
decrease. However, it will always increase on average 
and eventually the particles velocity will become much 
greater than the velocity of the wall. Then the work per 
cycle automatically becomes small and the conditions for 
the validity of the Fokker-Planck equation are satisfied. 

If the cavity is deformed while keeping its volume fixed 
then a very simple behavior emerges. In this case it 
has been shown [l6l - [l9j that the particles velocity dis- 
tribution becomes exponential irrespective of the con- 
tainer's shape and the deformation protocol, f{v, T)dv 



"/"^dv, where 



)t is proportional to time (num- 
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FIG. 2. An illustration of a single particle bouncing in a 
deforming cavity of constant volume. The driving protocol 
consists in repeatedly deforming the cavity between the two 
shapes shown. 

ber of cycles) and to the second moment of the velocity 
of the wall. Moreover if the cavity is sufficiently chaotic 
successive collisions are uncorrelated and the formalism 
holds even if the waiting time between cycles approaches 
zero (see for details). In this case the fundamen- 
tal equation for the energy distribution, P{E), assumes 
the Fokkcr-Planck form of Eq. (0) with A{E) = cE^I"^ 
and B{E) = cAE'^/'^ /{d + 1) where c contains informa- 
tion about the mass of the particle, the area or volume 
of the container and the velocity of the moving walls (see 
Rcf. |i3 for details). Using P{E) = [d - 2)/{2E) it is 
easy to verify that Eq. ([3|) is exactly satisfied with the 
correction included. Therefore, in the large d-limit, 
Eq. dl]) holds, and wc find a'^{E)/al^{E) — !• 2, consistent 
with the exponents a = 1 and s ~ 1/2. With the 
correction included we find: 

<^HE) ^ 2 + 3/d 

al^{E) l + l/d- ^> 

This result also follows from the exactly known single 
particle distribution [l^. This result was also derived 
for a single light impurity moving in a background of 
heavy atoms (Lorenz gas) . In this case a full microscopic 
description based on a Lorentz-Boltzmann equation is 
possible j20| . 

Weakly-interacting particles: Extending the above ex- 
ample, consider n weakly interacting particles in a de- 
forming cavity so that N = 2nd. We assume that during 
each cycle the particles can be treated as non-interacting, 
while between cycles the system rethermalizes at a fixed 
total energy, so that the velocity distribution of the par- 
ticles becomes Maxwell-Boltzmann rather than exponen- 
tial. A calculation similar to that of [l^, [l3l shows that 
A{E) and B{E) have the same functional form as in the 
non-interacting case but with different prefactors. More- 
over both A{E) and B{E) become extensive and the con- 
straint, Eq. ([3]), is satisfied with the 1/N correction be- 
coming negligible. Then Eq. ([1]) gives the asymptotic 
result cr"^ [E] / a1q{E) — > 2, consistent with the exponents 
a = 1 and s = 1/2. This result is identical to the sin- 
gle particle result in the large N limit, despite the very 
different single particle velocity distributions. 

In either noninteracting or interacting setups the two 
functions, A{E), B{E), can be experimentally obtained 



by measuring the first two cumulants of the work dis- 
tribution in one cycle: (w) = A{E) and {w^)c = B{E). 
Alternatively one can measure the average energy and its 
variance versus time and determine A{E) and B{E) as 
the slopes of these two functions respectively. 

Single particle in a time- dependent potential: Next we 
consider a classical particle in a harmonic trap, which is 
part of a larger system, e.g. a set of N identical particles, 
whose details define ^{E) and hence (i{E). We assume 
that the coupling to the rest of the system is weak and 
unimportant within the duration of a cycle, much like 
in the weakly-interacting particle gas example above. In 
contrast with the previous example, this setup illustrates 
driving a system with a local perturbation (which can, 
however, be applied independently to many different par- 
ticles). The particle's energy e between cycles is given by 



For simplicity we work in one dimension. For large N 
the probability distribution for (x, v) before the drive is 
p{x,v) cx exp (— /3(i?)e). Wc consider a driving process 
which consists of an impulse of magnitude F (x) At with 
At short enough so that the particle's position does not 
change appreciably during the drive. This assumption 
also guaranties that the coupling to the rest of the system 
is unimportant during a cycle. Under this impulse the 
velocity changes according to v ^ v-\-F {x) At. A, B are 
readily calculated and read 

A^(iF{x)f)At^, 

B-l{{F{x)f)At\ 

verifying fluctuation-dissipation relations ([3]). Taking 
/3 {E) oc E~°' and F {x) oc x^ we find 

A (X {x^'') (X e'' (X E""'. 

Using our previous convention, A E'^ ^ we see that s = 
ar. Different values of 77 = 2ar — 1 — a can be obtained 
using different impulse forces and systems. For example, 
for a system with a = 1/2, such as a Fermi liquid or a 
one-dimensional harmonic system, with r = 1 we obtain 
rj = —1/2 leading to the Gibbs-like regime with a'^/a'^^ = 
2. When r = 3/2 we are at the critical regime 77 = 0. 
Finally for r = 2 we obtain rj ~ 1/2 leading to the run- 
away regime where cr'^/cr'^q ^ E^/"^. 

II. CONCLUSIONS 

The main result of our paper is based on a fluctuation- 
dissipation relation connecting drift and diffusion of the 
energy in a driven system (Eq. ([3])). This relation is very 
closely connected to the recently discovered fluctuation 
theorems. In fact, in the supplementary material we give 
a rigorous derivation of Eq. ([3]) using the quantum version 
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of fluctuation theorems, which we extend to our setup of 
repeated cyche processes. In the original formulation the 
Jarzynski relation states that if a system starts from a 
Gibbs distribution the change in energy of the system, 
w, for cychc processes obeys the equality (e~^™) = 1. 
Here the angular brackets denote an average over both 
different realizations of the process and different initial 
conditions. As emphasized in Rcf. [l2| these relations 
holds very little information about the first few moments 
of the distribution of w, unless w is small. Only when 
it is small, a cumulant expansion of the Jarzynski rela- 
tion up to the second order in w recovers the fluctuation- 
dissipation relation ([3]) without the correction (see 
Supplementary Material for details). When the energy 
changes are large its moments are governed by details of 
the physical process which have to be accounted for (see 
p^). In this paper we (i) importantly, overcome the re- 
striction of small energy changes by considering a large 
change which is a results of many small changes (leading 
to Eq. ID)) and (ii) show that relation ^ is independent 
of the exact form of the initial distribution. 

We believe that some of the assumptions of the work 
can be further relaxed. In particular, it can be shown 
that relation (|3]) is valid for generic (non-cyclic) quasi- 
static process where the system is approximately in a 
steady state at each moment of time. In this case by A is 
related to the non-adiabatic part of the work w. Likewise 
it is plausible that the assumptions of complete relax- 
ation to a steady state between cycles are not necessary, 
at least for ergodic systems. Physically these assump- 
tions amount to a loss of correlations between different 
cycles which is inevitable in ergodic systems. The as- 
sumption of unitary dynamics can be also relaxed. For 
example, we can allow measurements of the energy dur- 
ing the protocol, which in quantum systems project the 
system to one of the energy eigenstates. Such measure- 
ments do not invalidate the derivation given in Appendix 
m These points will be addressed in a future study. 

The predictions of this paper can be experimentally 
tested in cold atom systems or in driven nuclear spins. 
For example, a trapping potential can be modulated to 
perform non-adiabatic work on the system. The result- 
ing energy distribution can then be probed using time of 
flight experiments. Likewise one can measure energy fluc- 
tuations due to a 1// noise in systems of neutral atoms 
near atom chips or in trapped ions in a setup similar to 
that discussed in Ref. 0. If there is no external cooling 
or dissipation in the system after its initial preparation 
(as is often the case) then the energy fluctuations induced 
by the noise should also agree with our predictions. 

Finally, we comment that energy fluctuations can be 
measured indirectly through averages or fluctuations of 
other observables, like the magnetization or correlation 
functions. In particular, it is easy to check that the ex- 
pectation value of every observable O and its variance, 
up to corrections, are given by 



(O) « O^c + 



"E '^mc 

2 dE^ ' 



60 
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mc 



[OeO 



2 2 

mc / 



where Omc is the microcanonical average of O at fixed 
energy E and SO^^ is the variance of O in the micro- 
canonical ensemble. Thus Eq. ([1]) will have implications 
to fluctuations of a wide class of observables in meso- 
scopic thermally isolated driven systems. 
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Boston University visitors program for its hospitality. 



III. METHODS 

We now outline a particularly simple derivation of 
Eqs. (H)) and (jS]). A more rigorous derivation based on 
fluctuation theorems and the unitarity of the evolution is 
given in the Supplementary Material. Our process con- 
sists of many repeated cycles during which the control 
parameter X{t) is varied in time, returning to its initial 
value at the end of each cycle. We assume that between 
the cycles the system reaches a steady state (or a diagonal 
ensemble [21[ in the quantum language) so that its state 
is fully characterized by its energy distribution. In er- 
godic systems this requirement can be satisfied by waiting 
between cycles a time which is longer than the relaxation 
time of the system. In non ergodic (integrable) systems 
this can be achieved by having a long fluctuating time be- 
tween cycles. This effectively leads to an additional time 
averaging which is equivalent to the assumption of start- 
ing from a diagonal ensemble. (For more details about 
relaxation to asymptotic states in integrable systems see 
Ref. and refs. therein). To make this discussion more 
concrete consider, for example, a compression and ex- 
pansion of the piston in Fig. [T] according to an arbitrary 
protocol. The gas is allowed to relax between the cycles 
(when the piston is stationary) at a fixed energy. For a 
weakly interacting ergodic gas such a relaxation implies 
that the momentum distribution of individual particles 
assumes a Maxwell-Boltzmann form together with a ran- 
domization of the coordinate distribution. For a nonin- 
teracting gas in a chaotic cavity the relaxation implies 
conservation of the individual energies of each particle 
and a randomization of the coordinates and directions of 
their motion. And finally for noninteracting particles in 
a regular non-chaotic cavity the relaxation implies a ran- 
domization of the coordinates within individual periodic 
trajectories. Therefore, in the beginning of each cycle 
there are no correlations between positions and veloci- 
ties of particles within the available phase space. 

If we make an additional assumption about ergodicity 
then the system between the cycles is fully described by 
the total energy. As we will see later, when we discuss 
specific examples, this assumption is not always neces- 
sary. Assuming that during each cycle a small amount 
of work is carried out on the system the energy distribu- 
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tion P{EA) can be described by the Fokker-Planck equa- 
tion ^ Q. The easiest way to derive the fluctuation- 
dissipation relation ([3]) is to note that under very gen- 
eral conditions the only attractor of Hamiltonian dynam- 
ics is a flat probability distribution for the occupation 
of different microstates (see Ref. [l^] for the classical 
case and Ref. [2^ for the quantum case), which is the 
maximum entropy state. Therefore the energy distribu- 
tion which is proportional to the many-particle density 
of states Ps{E) = CV,{E) should be stationary under 
the Fokker-Planck equation, implying that the current 
Js = -A{E)Ps{E) + ldE{B{E)Ps{E)) is a constant, 
which vanishes since Ps{E) = ioi E below the ground- 
state energy. Finally we use 13 (E) = Be Inri(ii^) to obtain 
Eq. ([3]). For a rigorous derivation of this result and the 
range of its applicability see Supplcmetary Material. 

The relation ([3]) allows us to make general statements 
about the energy distribution. In particular, the main 
result of the paper, Eq. ([IJ, immediately follows from 
Eqs. ^ and (jS]) to leading order in an expansion in 
To see this we first multiply Eq. ^ hy E and E'^ and 
integrate over all energies. In this way we obtain the 
differential equations describing the time evolution of {E) 
and (T^ = ^£;2^ _ {E)"^ , where angular brackets stand for 
averaging over P{E): 



dt{E) = {A{E)) 

dta^ = (S) + 2 i{A{E)E) 



{AiE)){E)). 



These two equations can be combined into a single one: 



da' {B)+2{{AE)~{A){E)) 



d{E) 



{A) 



(6) 



If the energy distribution P{E) is narrow, as in the case of 
large systems, we can evaluate the averages above using 
a saddle-point approximation and Eq. ([3]) . Then to order 



O(Af-i) we find: 
da' 



d{E) 



2r\{E)) 



(7) 



Integrating this equation immediately yields Eq. ([T]). 

Let us now comment on the regime of validity of 
Eq. ([T]). The derivation is based on the Fokker-Plank 
equation, Eq. ([2]), the generalized fluctuation-dissipation 
relation, Eq. and the saddle-point expansion in 

Eq. (O. The validity of the Fokker-Planck equation 
relics on the assumption that the work distribution is 
narrow and the average work for cycle is small. More 
specifically this equation is derived from a cumulant ex- 
pansion of the Crook's relation up to second order in 
the work (see Sec. U in Supplementary Material for de- 
tails). Let us here only mention two necessary con- 
ditions justifying the Fokker-Planck equation and the 
fluctuation-dissipation relation Eq. (i) The third 

(and higher order) cumulant of work per cycle are small 
l3'{E){w^)c < (w) = A{E). (ii) The average work per 
cycle is smaller than the product of temperature and the 
specific heat C„: l3{E){w) = 13{E)A{E) < C„. As ex- 
plained in the Supplementary Material (sec. U), if this 
condition is not satisfied there are corrections of order 
I3A'{E)/Cy to Eq. ©. Finally the saddle-point approxi- 
mation in Eq. ([7|) is justified if the energy fluctuations in 
the system are small. This is the case in large or meso- 
scopic extensive systems. 

We also note that our derivation implicitly relies on 
the assumption of ergodicity within the system. In par- 
ticular, we are assuming that P{E) is a differentiable 
function of energy. In integrable systems this is not nec- 
essarily the case [l5|. Then the validity of our results 



should be checked on a case by case basis. For exam- 
ple, as we showed above, Eq. ([3]) and the corresponding 
Fokker-Planck equation ([2]) describe the dynamics of a 
single particle in a cavity. In this case, however, the sec- 
ond term in the RHS of Eq. ([3|) is important and mod- 
ifles the width of the distribution even if we consider 
an ensemble of many nonintcracting (and therefore non- 
crgodic) particles. 
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SUPPLEMENTARY MATERIAL 

I. DERIVATION OF EQ. (4) FROM THE 
QUANTUM CROOK'S EQUALITY 

Here we will sketch the derivation of Eq. ^ relying 
only on the unitarity of the dynamics in the quantum 
case and the incompressibility of Hamiltonian dynamics 
(i.e. Liouville's theorem) in classical case. Our proof will 
be based on the Crooks theorem. As it was noted in pre- 
viously (see for example for a closed system obeying 
classical Hamiltonian dynamics, the Crooks theorem re- 
lies on the incompressibility of trajectories in phase space 
(Liouville's theorem) and microscopic time reversibility. 
Here we extend this proof to isolated quantum Hamilto- 
nian systems with a discrete spectrum. Our proof of the 
quantum Crook's relation bears some similarities with 
that discussed in Ref. , and is presented here for com- 
pleteness. This will emphasize some important proper- 
ties of the transition matrix, highlight that the Crook's 
relation does not rely on assumptions related to energy 
measurements, e.g. at intermediate steps, and extend 
the fluctuation-dissipation relation ([3]) to non-canonical 
distributions. 

Let us assume that a system prepared in a station- 
ary state, described by a diagonal density matrix in the 
energy basis, undergoes some process described by a uni- 
tary operator U (t) . According to standard quantum me- 
chanics the density matrix evolves in time according to 
p{t) = U{t)p{0)U^ (t). This means that the diagonal el- 
ements of the time evolved density matrix in the new 
energy basis arc given by 



T, 



(0) (8) 



where we used the fact that the initial density ma- 
trix is diagonal and introduced the transition probabil- 
ities T„i~in = |C^mnP (scc also Rcf. Q). The matrix 
Tm^n is doubly stochastic meaning that Tm^n = 
X]jn '^m^n = 1 • Wliilc thc first equality is simply the 
conservation of probability the second is a direct conse- 
quence of unitarity. It is easy to see that this equality is 
violated if there are losses in the system due to e.g. spon- 
taneous emission. Now let us imagine a time-reversed 
protocol described by the inverse evolution operator . 
From the definition of the transition probabilities it is 
clear that 



Tn- 



T„. 



(9) 



where T„_j.m refers to the reverse process. Let us 
comment that the transition probabilities also satisfy 
detailed-balance, Tn-^m = Tm^n, hi the two following 
situations: (i) if the Hamiltonian of the system is time- 
reversal invariant at each moment of time and the pro- 
tocol is time symmetric so that U{t) — U{T ~ t), where 
T is the period of the cycle, (ii) If the transition proba- 
bilities during one cycle are small and can be computed 



within first order in an adiabatic perturbation theory 0] , 
i.e. a perturbation theory in a basis evolving with the 
Hamiltonian (this theory also includes ordinary pertur- 
bation theory particular limit of small amplitude 
perturbations). Let us stress that detailed-balance only 
plays thc role in our proof for deriving the subleading 
OeB correction in thc relation ([3]). 

To proceed wc use thc energy distribution: 



P(S) =:^P„„,5(S-S„) 



(10) 



and relate transition probabilities between energy levels 
to the transition probabilities between energy shells: 

Te^E' = - En)5{E' - Era)T„^rn, (H) 

where Q{E) = ^„ S{E — En) is the many-body density of 
states. The factor l/il{E) ensures conservation of proba- 
bility: / (IE'Te^e' = 1- The master equation ([SJ is then 
given by 



P{E') ^ / dETE^E'PoiE) 



(12) 



We now multiply both sides of Eq. © by S{E — 
En)S{E' — Em) and sum over n, m to obtain 

n{E)TE^E' = n{E')fE'^E 

Denoting E' ~ E + w and using the fact that ^l{E) = 
cxp[S'(i?)] wc can rewrite the equation above as 

J-E~>E+we ' — 1e+w^E , (loj 

which is known as the Crooks relation (3]. 

To prove relation ([3]) we use Eq. (fTS]) . and expand the 
entropy and the transition probability Te^e+w in w- 



S{E + w)-S{E)kI3w- 



1 



W 



(14) 



^E — Te^E-w + wOeTe^E- 



as in the main text we assume a cyclic process. 

Note that when (w) is held constant and the system 
size is increased, the second terms on the RHS of each 
of the above equations scale as 1/A''. To leading order in 
1/iV, integrating Eq. (jTS]) over dw we obtain the Jarzyn- 
ski like relation <^e~''"') = 1, where the brackets represent 
an average over realizations of the process. Note that we 
did not make any assumption about an initial Gibbs dis- 
tribution. Taking the logarithm of this Jarzynski relation 
and performing a cumulant expansion we find 



[w 



from which Eq. ([3]) is obtained by using 

A={w), 



B 



[w 



(15) 



(16) 
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The condition for the vaUdity of this expansion is that 
the third cumulant of the work is small: 



(17) 



When the additional assumption of the detailed- 
balance holds we can use T = T in Eq. (fT3|) . We point 
again that detailed-balance is valid for arbitrary sym- 
metric protocols as well as for non-symmetric protocols, 
provided that the transition probabilities can be com- 
puted within first order of adiabatic perturbation the- 
ory. Since the work per cycle is assumed to be small, it 
is expected that the transition probabilities are also small 
and can be computed perturbatively. So the assumption 
o{ T = T between energy shells is likely to be generic. 
In particular, one can check that it is asymptotically sat- 
isfied at high energies for the piston example discussed 
in the main text even for asymmetric protocols. In this 
case, integrating Eq. p3p using the expansions we 
obtain 



exp 



2<, 



l-dE{w) 



(18) 



Taking the logarithm of both sides and performing the 
cumulant expansion of the exponent up to the order w'^ 
we find 



-dsiw). (19) 



It is easy to check that up to order 1/iV the equa- 
tions above imply Eq. ([3]): 2 A = f3B + dsB as long as 
{w)^/a^q is negligible compared to (3{w). Noting that 
(w^) — {'w'^)c + {w)'^ = B + this gives us a necessary 
condition of validity of relation 



A = (w) < Ta 



(20) 



Namely, the work per unit cycle should be small com- 
pared to the temperature multiplied by the specific heat. 
We note that even though wc derived Eq. ([3]) to the order 
of 1/N, it is actually correct to all orders in 1/A^. This 
relation is valid as long as the conditions p7|) and ([20)1 
are satisfied. 

Finally let us discuss extension of the fluctuation- 
dissipation relation ([3]) to arbitrary distributions. In or- 
der to do this we need to weight Eq. ([3]) with an energy 
distribution P{E) and integrate over energies. Then it 
is easy to check that in the Gaussian approximation we 
find 



(21) 



t In general situations, which can involve Berry phases, this state- 
ment is correct only in the appropriate adiabatic (co-moving) 
basis. 



This relation is clearly a generalization of Eq. (jSj. In 
particular, it reduces to Eq. ([3]) for the microcanonical 
distribution with — and it reduces to the conven- 
tional result obtained from the cumulant expansion of 
the Jarzynski relation (w) « /3/2{w'^)c for the canonical 
distribution a = acq- 



II. ADDITIONAL EXAMPLES 

In this section we will provide two additional exam- 
ples illustrating validity of our results to a classical and 
quantum interacting one-dimensional spin chains. 

A. An XY model in one dimension 

First we consider an XY-model on a one-dimensional 
lattice of size N. On each lattice site there is a sin- 
gle degree of freedom, which may be viewed as a two- 
dimensional unit vector. The interaction energy be- 
tween neighboring sites i,j is Hij = 1 — cos (Oij), where 
= Si — 6j is the difference between the angles at 
sites i and j. The total Hamiltonian is given by sum- 
ming over all nearest neighbors the interaction terms 
H = j) Hi.j- To drive the system, we assume that the 
angle at a specific site 9i is changed by a small amount 
59^ which is a fluctuating variable with zero mean, and 
other sites are unaffected. Such a protocol, for exam- 
ple, describes an interaction of the system with an ex- 
ternal fluctuating local magnetic field. The problem can 
be solved exactly in one-dimension. To do this, note 
that the change in energy of the system depends only on 
the two differences and Oi^i+i before the drive. The 

probability distribution of Oi i+i to order 0{1/N) is given 
by 



p(6lj,i+i) cx exp[/3(£') cos(6'^,i+i)]. 



(22) 



Using this expression, it is straightforward to calcu- 
late the values of the average work and its fluctuations 
A{E) = {w) ,B{E) = {w^)^ to order 



q2. 



A = {59'-) 1 - 



E 
N 



B 



2A 

T' 



(23) 



Here to relate E and /3 wc used the expression for the 
energy [3 E/N = 1 - [/i (/3) //q (P)], where /„ denotes a 
modified Bessel function of the first kind of order n. 

Substituting relations (j23p into Eq. ([T|) and numeri- 
cally integrating we obtain a^{E)/a'^(E) as a function 
of I3{E). The results are shown in Fig. [31 Note that there 
are two regimes. In the low energy regime A is to lowest 
order constant, s = 0, and a = 1, which gives 77 = —2, 
and cr^/cTei} — 1; ^ee Table HI The next order correction 
for small E/N can be obtained using Eq. ((^3)) and gives: 
a^/crlg ~ 1 — E/N. Then at high energies, close to infi- 
nite temperature, E/N = 1 and again /a'^g — >■ 1. This 
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FIG. 3. The ratio a^{E)/al^{E) as a function of E/N for a 
one-dimensional XY-model. The initial conditions are Eo — 
and ao = 0- 

is due to the finite phase-space available to this system, 
which allows the system to reach a stationary distribution 
with /3 = at a finite total energy density {E/N ~ 1). 

B. The transverse-field Ising model in one 
dimension 

Finally we consider a quantum transverse-field Ising 
model in one dimension described by the Hamiltonian: 

^=-E[5^." + ^I^I+i]+^'' (24) 
j 

where H' is a weak perturbation breaks integrability of 
the system but does not affect the dynamics during the 
driving protocol. For example, this perturbation can be a 
weak second nearest-neighbor spin-spin interaction. The 
assumption of smallness of H' is only important to make 
explicit analytic calculations. For simplicity we will con- 
sider the domain of non-negative values of the transverse 
field g. This system undergoes a quantum phase tran- 
sition at g = 1 Using a Jordan- Wigner transforma- 
tion the Hamiltonian assumes a quadratic form in the 
fcrmionic operators and can be fully diagonalized per- 
forming a Bogoliubov rotation in momentum space Q. 
The final Hamiltonian reads: 

^-T.4hh>^-ll (25) 

k 

where 

el = 2Vl + .g2-2,9Cos(fc) « 2^{1 -g)^ + k^ 

and 7^, arc quasi-particle creation and annihilation 
operators. Here in order to simplify analytic expressions 
we linearized the spectrum by taking g — I <^ 1 and the 
relevant momenta are much smaller than the ultra-violet 
cutoff given by the lattice: <^ tt. The ground state 
of this Hamiltonian, which is annihilated by all quasi- 
particle operators 7^., is factorized into momentum sec- 
tors. The excited states can be obtained by applying 



various combination of operators jj. to the ground state. 
If the external time-dependent perturbation is spatially 
uniform then due to momentum conservation only the 
excited states corresponding to pairs of excited quasi- 
particles with opposite momenta obtained by applying 
'Tfe'T-fe ground state participate in the dynamics. 

Moreover excitations to different momentum states are 
independent and the problem effectively splits into a col- 
lection of two levels-systems. Let us emphasize that this 
system does not become classical even in the infinite tem- 
perature limit. In disordered systems this feature was ex- 
plored e^. in the context of the many-body localization 
in Ref. [g. 

We now perturb the system by changing the ampli- 
tude of the transverse-field g{t) ~ gi + 5t{l — t/r), 
t G [0,t] from its initial value, gi, to an intermediate 
value, .92 = .91 + 5t/A and then back to gi. Here the pa- 
rameter 5 sets the velocity of the quench and 5t /A sets its 
amplitude. As a result of this process the occupation of 
the energy levels will change. Because different momen- 
tum modes are effectively independent from each other 
we will consider each two-level system separately. The 
presence of the weak integrability breaking perturbation 
ensures a Fermi-Dirac redistribution of the energy among 
different modes between different cycles. It is similar to 
the effect of a weak interaction between particles lead- 
ing to a Maxwell-Boltzmann single-particle distribution 
in the piston example or an assumption of weak cou- 
pling to the rest of the system in the single oscillator 
example discussed in the main text. The ergodicity of 
small weakly nonintegrable one-dimensional systems was 
recently tested numerically to a very good accuracy for 
hard-core bosons and fcrmions. which are closely related 
to the transverse field Ising model (see e.g. Ref. 0). 

Under these assumptions each cycle starts from the 
Fermi-Dirac distribution in each momentum mode. To 
avoid extra complications related to the singularities of 
the transition probabilities at the critical point we will 
additionally assume that dynamics occurs only in one 
phase, say g > \ (at finite temperatures we are inter- 
ested in, this assumption can be further relaxed). For 
slow quenches, if the rate 5 is small: 6 <C {gi — 1)^, such 
that the adiabaticity condition 5 ^ e^. is satisfied for all 
modes, the dynamics can be solved analytically. In par- 
ticular, in Ref. @ it was shown, that under these condi- 
tions the transition probability between the ground and 
excited states is approximately equal to (see Eqs. (20) 
and (87) in Ref. Hj) 

1 J^fc^ _ 2J^fc^ 

^""^ i2{k^ + {g,-lff- [el^r ^ > 

Note that there is an additional contribution to Eq. ([25)1 
which is rapidly oscillating at a frequency u ~ (32 — gi)T 
(see Eq. (20) in Ref. Q) and averages to zero either be- 
cause of adding contributions from different momentum 
modes or because of slight fluctuations of the quench time 
T from cycle to cycle. Actually the addition of the oscil- 
lating term to Eq. ([25]) will at most double the transition 



probability. As we will see below our results are not af- 
fected by the precise form of pk as long as pk is small, 
which is controlled by i5 and decays sufficiently fast with 
momentum k so that the work distribution remains suf- 
ficiently narrow. 

With these transition rates, the master equation for the 
occupation probabilities of different momentum states 
becomes particularly simple: 

Pi-Gr \ 1 - Pfc Pk \lpGr\ 

PkEx J \ Pk 1 - Pfe ; V PEcc J ' 

where p'^^ and refer to the probabilities of having 
a pair of fermions with momenta k and —k before and 
after the quench respectively and likewise, Pq^ and p^^, 
are the probabilities to have no fermions in the k and ~k 
momentum mode. Due to rethermalization of the system 
between driving protocols we have 

, _ exp[/?6f] _ cxp[-/?6f] 

2cosh[/3ef]'^^" 2cosh[/3ef] ^ ' 

Actually these expressions are somewhat modified due to 
presence of the inert modes with one fermion in either k 
or —k mode, but this modification does not affect our 
conclusions since the inert modes do not participate in 
transitions only effectively reducing the number of active 
modes by a factor of two in the regime of interest. From 
the master equation and the probabilities of occupying 
initial states above we find the average work and its sec- 
ond moment for individual momentum modes: 

(z«,)-26fp,tanh(/36f) 

{wl) = i2efrp, ^ 

From this equation we deduce the relation between 
A-k = (wk) and Bk ~ {w\)c'- 2Afe « jiBk is satisfied pro- 
vided that pfe ^ 1, which is the case for sufficiently slow 
quenches, and Pej} <^ 1 which means that the tempera- 
ture is big compared to the initial energy of the fermions. 
The latter condition is always satisfied if the tempera- 
ture is large compared to the gap in the system and that 
the relevant excited modes correspond to small momenta, 
which is the case for sufficiently slow quenches. If these 
conditions arc fulfilled then Eq. ([3]) is satisfied for the 
total work and its variance: 

k>0 k>0 

In Fig. |4] we show the calculated values of 2A and I3B 
for different initial temperatures for the protocol with 
gi — 1.1 and S = 0.05 (the time r drops out from the 
answer if St^ ^ 5i)- As the temperature increases the 
work per cycle. A, decreases due to fermion anti-bunching 
and the relation ([3]) is satisfied to a very good accuracy 
as soon as temperature becomes much bigger than the 
gap. 




FIG. 4. The validity of the Einstein-like relation for the 
transverse- field Ising model as a function of the initial temper- 
ature. The transition probability are calculated as in Eq. (|26p 
with 5 = 0.05 and gi = 1.1. In the inset we show the initial 
occupancy of energy level at T = 2.0. 



Let us make now a number of remarks on validity of 
our results and the importance of the protocol in quan- 
tum systems. The general analysis and derivation of the 
Fokker-Planck equation in the main paper as well as the 
relation in Sec. U of this Supplementary Material (see 
Ea. (|15p ) are based on a quantum formalism. Therefore 
as soon as the conditions f3{w) ^ and 0^{w^) <^ (w) 
are satisfied Eq. ([3]) should work. And indeed we illus- 
trated this here for a particular example of a transverse 
field Ising model. Instead of this model we could get e.g. 
fermions with two-body interactions and would come to 
similar conclusions. The smallness of work can be con- 
trolled by doing sufficiently slow quenches or by quenches 
of small amplitude similar to the classical case. There is 
an important subtlety related to the second condition of 
the smallness of the third and higher cumulants of work 
which distinguishes quantum and classical systems. In 
particular, in classical situations it is impossible to give 
a large energy to a system during a cycle unless the ex- 
ternal parameter changes suddenly. E.g. in a deform- 
ing cavity example during each collision with a wall the 
particle can gain at most velocity 2V, where V is the 
velocity of the wall. This means that in classical situa- 
tions the work distribution is typically bounded and the 
smallness of average work usually implies smallness of 
its cumulants. In quantum systems the situation is very 
different. Namely in any protocol it is possible to give 
the system arbitrarily large energy. For smooth proto- 
cols, where the external parameter changes analytically 
in time the transition probability to high energy states 
decreases exponentially like in the conventional Landau- 
Zener problem and such transitions do not affect cumu- 
lants of work. However, for non-analytic protocols where 
e.g. amplitude, velocity (like in our case) or accelera- 
tion experience a discontinuity the transition probability 
to high energy states decreases only algebraically with 
energy like in Eq. (|26p (see Ref. for more details). 



This means that high enough cumulants of work neces- 
sarily diverge. Of course real experimental protocols are 
always smooth and these divergencies are cutoff, how- 
ever the degree of smoothness introduces a new quantum 
scale into the dynamics. Thus the work distribution can 
be effectively wide or narrow in the sense of satisfying 
the condition f3'^{w^) ^ (w) depending on the ratio of 
this new energy scale and temperature. 

Since this discussion is not directly related to this pa- 
per we postpone it until future work. Here we only explic- 
itly analyze another protocol where the coupling changes 
quadratically in time: g{t) = Qi + {e/2)t'^{l — t/r)'^. In 
this case in the slow limit e ^ (ffi — l)"^ instead of Eq. 
we get (see Ref. [1] for details): 



~ 2 (4^- ^^^^ 



As we see indeed the transitions to the higher energy 
states are suppressed even more than for the linear pro- 
tocol so that more cumulants of work now converge at 
small energies. As a result the relation ([3]) is satisfied 
even at smaller temperatures, see Fig. [S] In Fig. [5] wc 
show the values of 2A and f3B obtained using Eq. ((29)) 
with gi = 1.1 and e = 0.024 (the time r drops out from 
the answer if er^ ^ gi). With this choice of the pa- 
rameters the peak value of pk in the two protocols (lin- 
ear, Eq. (p6)) . and quadratic, Eq. ((29)) ) are equal. The 
Einstein-relation, 2 A = (3B, is satisfied when the relative 
error {B/T - 2A)/2A becomes less than 5%. With this 
definition the Einstein-relation is satisfied at T > 2 for 
the linear and T > 0.8 for the quadratic protocol. 




FIG. 5. Einstein-like relation for the transverse-field Ising 
model as a function of the initial temperature. The transition 
probability are calculated as in Eq. (|29|) with e = 0.024 and 
gi = 1.1. In the inset we show the initial occupancy of energy 
level at T = 0.8. 

In the regime of validity of the fluctuation-dissipation 
relation ((3|) we can expand expression ((28p to first order 
in /3 to obtain A ^ /3 ^ for both the linear and 
quadratic protocols described above. Moreover, if the 
temperature T is in turn much smaller than the cutoff en- 
ergy scale given by J (which is nothing but the Fermi en- 
ergy in the fermion representation) then E ^ , which is 
the case for weakly interacting fermions if T ^ Ep- From 
these considerations it immediately follows that a = 1/2, 
s = -1/2 and ry = -5/2. Therefore a^/crlq 2/5. At 
high energies, i.e. close to the infinite temperature limit, 
A — and /a1^ must tend to 1. As in the XY-model 
this is due to the finite size of the Hilbcrt space (see 
discussion after Eq. ([23|))- This example shows that for 
weakly interacting fermions (to which the Ising spin chain 
is equivalent) it is easy to get distributions significantly 
narrower than canonical without any special fine tuning. 
We expect this to be generic result following from the 
Pauli exclusion principle and suppression of transitions 
at higher temperatures (see Ref. [1] for additional dis- 
cussion). 
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